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Abstract 

Based on the work of Green, Tao and Ziegler, we give asymptotics when 
N — > oo for the number oftixn magic squares with their entries being prime 
numbers in [0, N]. For every n > 3 we give appropriate systems of linear forms 
(or equivalently basis) describing all n x n magic squares with integer entries 
and we calculate the complexity of these systems in the Green and Tao sense. 
We compute the precise asymptotics for the cases n — 3 (complexity 3) and 
n = 4 (complexity 1), and the given algorithm works for n > 5 (complexity 1). 
Finally, we show that the asymptotics are exactly the same if we impose that 
all the entries of the magic squares have to be different. 



1 Introduction 

In this article we will be interested in the asymptotic number of magic squares 
with their entries being prime numbers in [0, N] when N — > oo. 

The generalization of an important conjecture by Hardy and Littlewood, says 
(Conjecture 1.4 in [5]): 

Conjecture 1. (Generalised Hardy-Littlewood conjecture) Let N , d, t, L be positive 
integers, and let Vl/ : Z d — > Z' ; ^ = (ipi, . . . ,tpt)j be a system of affine-linear forms 
(none of the ipi is constant and no two of them are rational multiples of each other) 
with size \\^\\n ■= Yn=i Yfj=i l^i( e 3')l+Z)i=i - L > where eachipi = ^+^(0) 

and e±, . . . , ed is the standard basis for Z d . Let K C [-N, N] d be a convex body. Then 
we have: 

l*nz<n •-'(*•)! = (1 + 0,^(1))^ n * + °<.«(,5v)- 

b p prime \ o / 

where P = {2,3,5,. . .} denotes the prime numbers, the archimedean factor (3^ is 
vold(K n^!' 1 (M} + )) and the local factors f3 p are defined as f3 p = K ngZ d Yl i Az p {ipi(n)) 



for each prime p with A% p being the local von Mangoldt function, A% p : Z p — > M + , 
A Zp (0) = and A Zp (b) = jfc = ^ if b ± 0. 

Informally speaking, in the same sense that the Prime Number Theorem states 
that if a random integer is selected near N its probability of being prime is asymp- 
totically 1/logiV, the conjecture asserts that the probability that a randomly se- 
lected point in ^ f (Z rf )nZ^_ "of magnitude N" has prime entries in all its coordinates 
is asymptotically Y[ p /W log* N. 

Green and Tao introduced in jS] a "measure of how complicated such a problem 
is". It is called complexity (or Cauchy-Schwarz complexity): 

Definition 1.1. If = (tpi, . . . ,ip t ) is a system of affine-linear forms, it has i- 
complexity at most s if one can write the set oft—1 forms {tp\, . . . , ipi-i, tpi+i, • • • , ipt} 
as a union of s + 1 sets, in such a way that tpi does not lie in the affine-linear span 
over Q of any of these sets. The complexity of \£ is the least s for which the system 
has i-complexity at most s for all 1 < i < t, or oo if no such s exists. 

The Main Theorem of Green and Tao in that same article, says: 

Theorem 2. Suppose that the Inverse Gowers-norm conjecture GI{s) and the 
Mobius and Nilsequences conjecture MN(s) are true for some finite s > 1. Then 
the generalised Hardy-Littlewood conjecture is true for all systems of affine-linear 
forms of complexity at most s. 

In [9] the Mobius and Nilsequences conjecture MN(s) is proved for every s and 
in |10j the Inverse Gowers-norm conjecture GI(s) is proved for every s. Then, 
Conjecture [T] is in fact a theorem for every system of affine-linear forms of finite 
complexity. 

For the problem of n x n magic squares, as we shall soon see, we will be dealing 
with linear instead of affine-linear forms. In Section [2j we will find suitable linear 
forms for our problem for every n > 3. In Section [3j we will calculate the complexity 
of all these systems of linear forms. In Section [4j we will deal with the computation 
of the volumes of the polytopes that arise in our problem. In Section [5j we will face 
the computation of the local factors (3 p . Based on the work of the other sections, 
in Section [6] we will be able to give our asymptotics. Finally, in Section [7j we will 
show that the asymptotics are the same if we impose the condition that our magic 
squares must have different entries. 
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2 Basis for Magic Squares 



Definition 2.1. For any n > 3, an rex n M-magic square is an nx n array with its 
n 2 entries being real numbers, in such a way that the sum of the elements in every 
row, column or any of the two main diagonals is the same, the magic sum S. If all 
its entries are integers, we will say it is a Z-magic square. 

Remark 2.1. For every n > 3, n x n M.-magic squares have the structure of a M.- 
vector space and rex n Z-magic squares form an abelian group with the sum, and so 
a Z-module, and we will soon see that we can find a basis for this Z-module. Most 
of the time in this article we will be only interested in Z-magic squares. When we 
want to specify that the re 2 entries of a magic square are different ( in fact, what is 
usually called a magic square is a Z-magic square with different entries), we will 
explicitly state it. Also, we will specify when we want the entries to belong to certain 
set (e.g. {0, . . . ,N} or the primes in this set). 

For a 3 x 3 Z-magic square, naming its entries x\,X2, ■ ■ ■ ,xg from left to right 
and from top to bottom (as we will always do in the rest of the articl^]), we have: 

(xi + x 5 + xg) + (x 2 + x 5 + x 8 ) + (x 3 + x 5 + x 7 ) = 35 

and since x\ + x 2 + £3 = xj + x$ + xg = S, we have that S = 3x^. So, if the integer 
£5 = a, then x\ and xg will be respectively + 6 and a — b for some integer b. In 
the same way, X3 and X7 will be respectively a + c and a — c for some integer c. The 
other four entries of the square are now determined, as we can see in Figure [T] 

xi X2 X3 a + b a — b — c a + c 

xa 25 ^6 = a — b + c a a + b — c 
X7 x$ xg a — c a + b + c a — b 

Figure 1: Entries of a 3 x 3 magic square. 

So, we already have our system of linear forms for 3x3 Z-magic squares, : 
1? — > Z 9 , VP(a, 6, c) = (a + b,a — b — c,a + c,a — b + c,a,a + b — c,a — c,a + b + c,a — b). 

Not only for any a, b, c £ Z we have that ^(a, b, c) gives the entries of a Z-magic 
square, but also, because of the way we constructed it, we have that every 3x3 

1 We will think about these vectors with n 2 coordinates both as column vectors and as n x n 
squares. So, we will sometines make an abuse of notation refering to the second row of a vector 
(entries n+ 1, . . . , 2n) or its top-left to bottom-right diagonal (entries 1, n + 2, In + 3, . . . , n 2 ). 
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Figure 2: Z-basis for 3x3 Z-magic squares. 

Z-magic square is obtained for some a,b,c G Z (and this last thing is important 
because, if not, our linear forms could be describing only a subset of all Z-magic 
squares). In other words, and this approach will be useful for magic squares with 
larger side, looking at the linear forms, we conclude that given any 3x3 Z-magic 
square, we can obtain it as an integer linear combination of the three magic squares 
given in Figure [2j 

Then, we have a basis for the Z-module of 3 x 3 Z-magic squares. Obtaining this 
basis is clearly equivalent to obtaining the linear forms we are interested in: if we 
write the entries of each one of the squares of the basis as column vectors, one after 
another, and we look at the rows, we will read the coefficients of the linear forms. 

Since a finite set of vectors with integer entries is linearly independent over Z if 
and only if it is linearly independent over M. (the "if" part is obvious and the "only 
if" part can be deduced from the fact that, thinking the vectors of the set as rows 
of a matrix, elementary row operations over Z give a row reduced form over Z of 
the matrix), in what follows we can interpret if we want "linearly independent" as 
"linearly independent over W . 

Definition 2.2. We will say that a linearly independent set ofnxn Z-magic squares 
is a Z-basis for nxn Z-magic squares if every nxn Z-magic square can be obtained as 
an integer linear combination of its elements. We will say that a linearly independent 
set of nxn M-magic squares is an M.-basis for nxn M-magic squares if every nxn 
M-magic square can be obtained as a real linear combination of its elements. 

We remark that the given Z-basis for 3x3 magic squares (as every other Z-basis) 
is also an M-basis. But not every M-basis (even if all its elements have only integer 
entries) is a Z-basis. So, how to obtain an M-basis for nxn magic squares and how 
to make sure it is a Z-basis? 

First of all we can easily obtain a system of linear equations that defines nxn 
M-magic or Z-magic squares. First, we have n — 1 conditions (linear equations) of 
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Figure 3: Matrix A for n = 4. 



type A saying that the sum of the entries of the i-th row is equal to the sum of the 
entries of the z + 1-th row (for i = 1, . . . , n — 1). Then, we have other n — 1 conditions 
of type B saying that the sum of the entries of the i-th column is equal to the sum 
of the entries of the i + 1-th column (for i = 1, . . . , n — 1). These conditions also 
imply that the sum of every row is equal to the sum of every column. Now, we have 
condition C saying that the sum of the entries of one main diagonal is equal to the 
sum of the entries of the other main diagonal. Finally, we have condition D saying 
that the sum of the entries of the main diagonal from top-left to bottom-right is 
equal to the sum of the entries of (say) the first column. These 2n conditions clearly 
define n x n M-magic or Z- magic squares. If we name x the column vector with 
entries x\, . . . ,x n 2, and the column vector with 2n zeros, then n x n R- magic or 
Z-magic squares are defined by the system of linear equations Ax = 0, where A is 
the matrix formed with the coefficients of the 2n conditions previously described, 
which, for example, for n = 4 would be the one given in Figure [3} 

Are the 2n rows of the matrix A linearly independent? Yes, they are. An 
elegant way to see it is to exhibit, for every i = 1, . . . , 2n, a vector that satisfies the 
conditions in all the previous rows but not the condition in row i. The vectors of 
type a with ones in the first in positions (the first i rows) and zeros in the rest (for 
i = 1, . . . , n — 1) satisfy all type A conditions except from condition i (in particular 
all conditions above row i and not condition in row i). The same happens with 
vectors of type b that have ones in the positions "j + multiple of n" for j = 1, . . . , i 
(the first i columns) and zeros in the rest (for i = 1, ... ,n — 1): they satisfy all 
conditions above n — 1 + i but not condition n — 1 + i. The vector c with one in 
all its positions except from positions 1, n + 2, 2n + 3, . . . , n 2 (the main diagonal 
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Figure 4: Examples of vectors of type a and type b and vectors c and d for n = 5. 

from top-left to bottom-right) where it has zeros, satisfies all conditions of type A 
and type B but not condition C. Finally, the vector d with zero in all its positions 
except from positions 1, n, n + 2, 2n — 1, 2n — 3, 3n — 2, . . . , n 2 — n + 1, n 2 (the two 
main diagonals) where it has ones (or a two in position in 2 + l)/2 when n is odd), 
satisfies all conditions of type A and type B and condition C but not condition D. 
Figure [4] shows an example of how vectors of type a and type b and vectors c and d 
look like when we write them asnxn squares, for n = 5. 

We deduce that the rank of matrix A is 2n, and so the R- vector space of R- magic 
squares has dimension n 2 — 2n. To find an R-basis of this vector space we need to 
find n 2 — 2n linearly independent vectors that satisfy the In equations, that is, 2n 
linearly independent magic squares. This can be done easily for every fixed n, since 
we know how to solve a system of linear equations. But we also want our R-basis 
to be a Z-basis. In principle it would not be obvious even that a Z-basis exists. 

We could, for example, look for an appropriate matrix obtained from A by 
elementary row operations, with 2n of its columns being the vectors of the canonical 
basis of the 2n-dimensional Euclidean space and with all the rest of its entries being 
integer numbers. Instead of that, and maybe more naturally, we will give an explicit 
Z-basis for n x n magic squares and this will prove in particular that such a matrix 
exists for every n > 4. 

We know the number of linearly independent vectors that we need: n 2 — 2n. In 
the case n = 3 we had 3 2 — 2 • 3 = 3 vectors in our Z-basis. For n = 4 we need 
4 2 — 2 • 4 = 8 linearly independent magic squares that form a Z-basis. 

The easiest way to understand our basis for n x n Z-magic squares for n > 4 
is the next. We will carefully select some positions of the square, more concretely 
n 2 — 2n positions, in such a way that one time we fix the entries in these positions, 
the rest of the entries of the magic square are determined (we have discovered that 
this idea is already pointed out in [15J). The selected positions will form the skeleton 
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of the basis (following the notation in |15j). Then, the squares of the basis will be 
the only n 2 — 2n squares that have one as an entry in one position of the skeleton 
and zeros in all the rest of it. 

For n = 4, selecting the next skeleton of 8 positions, 



we have the Z-basis for 4x4 Z-magic squares shown in Figure [5j 

1000 0100 0010 

0001 0001 0001 

0100 0010 0010 

0010 1000 11 -1 
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Figure 5: Z-basis for 4x4 Z-magic squares. 



So everything has worked well but, in principle, we have been lucky. There 
is not an obvious reason why the entries outside the skeleton should be integers. 
As an example showing that this does not always happen, a 3 x 3 IR-magic square 
is characterised by the three entries of its first row, and the only M-magic square 
starting with one, zero and zero in these three positions is: 

1 

-2/3 1/3 4/3 
2/3 2/3 -1/3 . 



With this in mind, we define our basis for n > 5: 
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Definition 2.3. (Elephant basis) For n > 5, consider the skeleton in an n x n 
square consisting of all the elements of the first row, all the elements of rows 2 to 
n — 2 except from the last one of each row and all the elements of row n — 1 except 
from the second, the (n — 2)-th and the last one. The nx n elephant basis is formed 
by the n 2 — 2n nxn Z-magic squares that have one as an entry in one of the positions 
of the skeleton and zeros in all its other positions. 



Figure 6: Skeleton of the elephant basis. 

To see that the definition is right and that it gives us what we want, we prove 
the next lemma. 

Lemma 3. For every n > 5 the elephant basis is a Z-basis for nxn TL-magic 
squares. 

Proof. First of all we will show that if we fix any integer entries in the positions of 
the skeleton (in particular if one of them is 1 and all the rest are 0) , they determine 
a unique M-magic square which is also a Z-magic square. We first observe that the 
skeleton contains the first row of the square, so adding the values in these positions 
gives the magic sum of the square. Since the magic sum is an integer number, we can 
now determine the integer numbers of the last position of rows 2 to n — 2 (because of 
the sum on their rows), the last position of columns 1, 3 to n — 3 and n — 1 (because 
of the sum on their columns) and the one in the bottom-right corner of the square 
(because of the sum on its main diagonal). With these integers already determined, 
we can determine the integer in position 2 of row n — 1 (because of the sum on its 
main diagonal) and the one in the last position of the same row (because of the sum 
on its column). Now, we can determine the integer in the second position of the last 
row (because of the sum on its column) and the one in position n — 2 of row n — 1 
(because of the sum of its row) . The value in position n — 2 of the last row is the 
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only one left to be determined. Since the sum of all the elements of the square is 
the same if we sum them by rows or by columns, we will of course obtain the same 
integer if we determine it with the values of column n — 2 or the ones in the last 
row. 

Now, the vectors are linearly independent by construction (each of them has a 
1 in a position where all the others have a 0), and, given a Z-magic square, we can 
look at the values it has in the positions of the skeleton and we can express it as 
the linear combination with these coefficients in the corresponding vectors of the 
elephant basis. So the elephant basis is a Z-basis. □ 

Following exactly the procedure described in the proof of Lemma [3] we can de- 
termine all the entries of all the squares in our elephant basis. If we write them 
as column vectors, one after another, forming a matrix, the rows will give us the 
linear forms we were looking for. Figure [7] gives this n 2 x (n 2 — 2n) matrix. We will 
identify each row of coefficients with the corresponding linear form. The n 2 — 2n 
forms marked in gray in the figure are the "canonical" linear forms with all their 
coefficients equal to zero except from one which is 1 and correspond with the posi- 
tions of the skeleton of the elephant basis. We will call these the trivial linear forms. 
The other 2n forms will be called the nontrivial linear forms. 

Of course, the first n vectors (columns) have magic sum 1 and all the others 
have magic sum 0. The pattern on each row, if we look at the blocks of n — 1 
columns indicated (from column n + 1 to column n 2 — 3ra + 3), should be obvious. 
Maybe only some words about the patterns on the third and (n + 3)-rd row from the 
bottom should be said. We explain what happens in the third row from the bottom 
(the other is very similar) . The first n and the last n — 3 positions do not give any 
problem, so we look at the n — 3 blocks of n — 1 positions from column n + 1 to 
column n 2 — 3n + 3. To obtain the second block from the first (observe that in the 
case n = 5 the first block is 2 2 0), we substract 1 in the second position and add 
1 in the third and also substract 1 in the (n — 2)-nd and add 1 in the (n — l)-st. 
To obtain the third block from the second, the —1 and +1 on the left move one 
position to the right (we substract 1 in the third position and add 1 in the fourth) 
and the —1 and +1 on the right move one position to the left (we substract 1 in 
the (n — 3)-rd position and add 1 in the (n — 2)-nd). And we continue in the same 
fashion. 
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As we announced in our comment on page |6j obtaining n 2 — 2n trivial linear 
forms proves in particular that there is a matrix obtained from A by elementary 
row operations, with 2n of its columns being the vectors of the canonical basis of 
the 2n-dimensional space and with all the rest of its entries being integer numbers. 
In fact, we can very easily write such a matrix explicitly from the linear forms: the 
vectors of the canonical basis of the 2ra-dimensional space will be in the In columns 
of the matrix corresponding to the nontrivial linear forms and the other n 2 — 2n 
entries of the rows of A will be the opposites of the coefficients of the nontrivial 
linear forms. 

Then we have our linear forms for every n. Now what? We turn to the complexity 
problem. 

3 Complexity 

Althought the results mentioned in Section [T] tell us that we can use Conjecture 
[I] as a theorem for every system of affine-linear forms of finite complexity, deter- 
mining the exact complexity of a problem is theoretically important. Problems of 
complexity 1 were essentially present in the work of Hardy-Littlewood and Vino- 
gradov, though not in this language. To solve problems of complexity 2 we need to 
use the results in [6], [7] and [8]. To solve problems of complexity 3 or higher we 
need to use the results in jS], [S] and [TU] . 

We have Z-basis (or equivalently linear forms) for 3 x 3, 4 x 4 and n x n (for 
n > 5) Z-magic squares in Figure [2j Figure [5] and Figure [7] respectively. We will 
calculate the complexity (see Definition of these systems of linear forms. 

3.1 Complexity for 3x3 Z-magic squares 

The system of linear forms given by the Z-basis of Figure [2] is, as we pointed 
out, \P(a, 6, e) = (a + 6, a — b — c,a + c,a — b + c,a,a + b — c,a — c, a + b + c, a — b). We 
will indentify these linear forms with the row vectors (1, 1, 0), (1, —1, —1), (1, 0, 1), 
(1,-1,1), (1,0,0), (1,1,-1), (1,0,-1), (1,1,1) and (1,-1,0). 

Lemma 4. The complexity of the system of linear forms given by the Z-basis for 
3x3 7L-magic squares of Figure [#] is 3. 

Proof. First of all, we show that the i-complexity of every form is at most 3. In 
order to see it, given the form i, we split the others into four sets in such a way that 
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the form i cannot be written as a linear combination over Q of the forms in any of 
the four sets. We write explicitly the four sets for every form: 

• (l,l,0)-complexity at most 3: {(1, 0, 0), (1, 0, 1)}, {(1, -1, 0), (1, 0, -1)}, 
{(1,1,1), (1,-1,1)}, {(1,1,-1), (1,-1,-1)}. 

• (l,-l,-l)-complexity at most 3: {(1, 1, 0), (1, -1, 0)}, {(1, 0, 1), (1, 0, -1)}, 
{(1,0,0), (1,-1,1)}, {(1,1,-1), (1,1,1)}. 

• (1, 0, l)-complexity at most 3: {(1, 0, 0), (1, 1, 0)}, {(1, -1, 0), (1, 0, -1)}, 
{(1,1,1), (1,-1,-1)}, {(1,-1,1), (1,1,-1)}. 

• (1, -1, l)-complexity at most 3: {(1, 1, 0), (1, -1, 0)}, {(1, 0, 1), (1, 0, -1)}, 
{(1,1,1), (1,1,-1)}, {(1,-1,-1), (1,0,0)}. 

• (l,0,0)-complexity at most 3: {(1, 1, 0), (1, 0, 1)}, {(1, -1, 0), (1, 0, -1)}, 
{(1,1,1), (1,-1,1)}, {(1,1,-1), (1,-1,-1)}. 

• (1, 1, -l)-complexity at most 3: {(1, 1, 0), (1, -1, 0)}, {(1, 0, 1), (1, 0, -1)}, 
{(1,1,1), (1,-1,1)}, {(1,-1,-1), (1,0,0)}. 

• (1, 0, -l)-complexity at most 3: {(1, 0, 0), (1, 1, 0)}, {(1, -1, 0), (1, 0, 1)}, 
{(1,1,1), (1,-1,-1)}, {(1,-1,1), (1,1,-1)}. 

• (1, 1, l)-complexity at most 3: {(1, 1, 0), (1, -1, 0)}, {(1, 0, 1), (1, 0, -1)}, 
{(1,0,0), (1,-1,1)}, {(1,1,-1), (1,-1,-1)}. 

• (1, -1, 0)-complexity at most 3: {(1, 0, 0), (1, 0, 1)}, {(1, 1, 0), (1, 0, -1)}, 
{(1,1,1), (1,-1,1)}, {(1,1,-1), (1,-1,-1)}. 

So the complexity of the system is at most 3. Now we will prove that, for 
example, the (1, 0, 0)-complexity has to be 3 and no less and this will prove that the 
complexity of the system is 3. We suppose for a contradiction that we can split the 
other eight forms in three sets in such a way that (1, 0, 0) is not in the linear span of 
any of these sets. Since (1, 0, 0) is in the linear span over Q of (1, 0, 1) and (1, 0, —1), 
these two forms have to be in a different set, and the same happens with (1,1,0) 
and (1, —1,0). Since we only have three sets, two of these four forms have to be in 
the same set, so we have that one form of the shape (1, *, 0) and one of the shape 
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(1,0,*) are in the same set, where * can be 1 or —1. In that set we cannot have 
any other form, since (1, 0, 0) would be in the linear span of them together with any 
(1, *, *). So the first of our three sets is of the form {(1, *, 0), (1, 0, *)}. Now, since 
(1, 0, 0) is in the linear span over Q of (1, 1, 1) and (1, —1, —1), these two forms have 
to be in a different set, and the same happens with (1, —1, 1) and (1, 1, —1). Then, 
the two forms of the shape (1, *, 0) and (1, 0, *) that are not in our first set, must be 
in different sets, since if they are in the same set they would be the only two forms 
in that set and then we would have more than three sets. Then a form of the shape 
(1,*,0) is in the second set and a form of the shape (1,0,*) is in the third. Now, 
either (1, 1, 1) and (1, —1, 1) are in one of these sets and (1, —1, —1) and (1, 1, —1) in 
the other or (1, 1, 1) and (1, 1, —1) are in one these sets and (1, —1, —1) and (1, —1, 1) 
in the other. In the first case, (0, 1, 0) is in the linear span of both sets, and (1, 0, 0) 
in the linear span of the set with the form (1, *, 0), which is a contradiction. In the 
second case, (0, 0, 1) is in the linear span of both sets, and (1, 0, 0) in the linear span 
of the set with the form (1,0, *), which is also a contradiction. □ 

This tell us two things, that complexity can be somewhat complicated and that 
for 3x3 squares we need very "strong" and recent results. 

3.2 Complexity for 4x4 Z-magic squares 

We write the vectors of the Z-basis of Figure [5] as columns of a matrix (see 
Figure [8]) and, looking at the rows, we have the coefficients of the 16 linear forms, 
ipi, ... , ipiQ, we are interested in. We identify each form with its corresponding row 
vector in 1? . Because of the way we constructed the basis, observe that ipi, . . . ,ip7 
and ijjg are the vectors of the canonical basis of the 8-dimensional space. 

We analyse the ^-complexity of each form: 

• ipi, tps, V^io and V'is-complexities are at most 1 because ipi, ips, V^o and ipis 
are linearly independent and not in the linear span of the rest (all of them 
have a value different from in their first position and all the rest have a zero 
in that position) . 

• ip2, ip3, ipn, ipi3 and ^^-complexities are at most 1 because tp2, i>3, ips, tpii, 
■013) V*14 and -015 are linearly independent and not in the linear span of the 
rest (all of them have a value different from in their second or third positions 
and all the rest have a zero in both of these positions) . 



13 



/ 1 























\ 


(1 


1 




























1 




























1 




























1 




























1 




























1 







1 


1 


1 


1 


-1 


-1 


-1 




























1 




1 








-1 


1 





1 


1 







1 


1 


2 


-1 


1 





-1 



















1 


1 


-1 







1 


1 


1 


-1 








-1 










1 


2 


-1 


1 


1 


-1 




1 





-1 


-1 


1 


1 


-1 


1 




\ o 








-1 


1 








1 


/ 



Figure 8: Matrix with the vectors of the elephant basis as columns for n = 4. 

• ^4 and ^i6-complexities are at most 1 because ifii, ipg, tpio, ipu, ipis, ipu, tpi^ 
and tpiQ are linearly independent and not in the linear span of the rest (all 
of them have a value different from in their fourth position and all the rest 
have a zero in that position). 

• ips and ^g-complexities are at most 1 because ips, ipQ, ipj, ipg, -011; ipi3, V'u 
and V'ls are linearly independent and ips and ip§ are not in the linear span of 
the other eight forms. 

• ipQ and ^^-complexities are at most 1 because ip^, V>8> ^llj "012) 4>ia and ■(/'is 
are linearly independent and not in the linear span of the rest (all of them 
have a value different from in their sixth position and all the rest have a zero 
in that position). 

• V^-cottipl^xity is at most 1 because ip?, ipg, tpio, ip\2, ?pi4 and 7/>i5 are linearly 
independent and not in the linear span of the rest (all of them have a value 
different from in their seventh position and all the rest have a zero in that 
position). 

So the ^-complexity of every form is at most 1 (and it cannot be less since the 
16 linear forms are not, of course, linearly independent). Then: 

Lemma 5. The complexity of the system of linear forms given by the TL-basis for 
4x4 "L-magic squares of Figure^ is 1. 
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3.3 Complexity for n x n Z-magic squares for n > 5 

Recall that the "white" rows in Figure [7] give the coefficients of the nontrivial 
linear forms of the elephant basis for every n > 5. 

Lemma 6. By elementary row operations in the 2n x (n 2 — 2n) matrix with the 
coefficients of the nontrivial linear forms of the elephant basis as rows (matrix formed 
by the "white" rows in Figure^!, the matrix of Figure^ can be obtained. 

Proof. The rows of the matrix 
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list the coefficients of the linear combination of the rows in the matrix of coefficients 
of the nontrivial linear forms that give the rows of the matrix in Figure [9j □ 

Proposition 7. The complexity of the system of linear forms that define n x n 
Z-magic squares given by the elephant basis for n > 5 is 1. 

Proof. From the previous result, and looking at Figure [9| we know that the 2n 
nontrivial linear forms of the elephant basis are linearly independent (the 2n columns 
marked in gray have zeros in every position except from one). 

Given any trivial linear form different from the ra-th one, it is impossible to write 
it as a linear combination of the rows of the matrix in Figure [9] (and so as a linear 
combination of the nontrivial forms): If such a linear combination exists and we 
look at the positions corresponding to the columns in grey in Figure [9] at most one 
of them will be different from zero. This implies that all the coefficients of the linear 
combination except from at most one are zero. Now, if the trivial linear form is 
different from the n-th one, it cannot be obtained as a linear combination with only 
one nonzero coefficient of the rows in the matrix of Figure [9] because none of them 
is proportional to a trivial linear form. 
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If we join the first and third trivial linear forms ((1, 0, . . . , 0) and (0, 0, 1, 0, . . . , 0)) 
to the nontrivial linear forms, we have a set of 2n+2 linearly independent forms: The 
trivial form (0, 0, 1, 0, . . . , 0), say, is linearly independent with all the nontrivial ones 
by the last paragraph argument. Now, suppose for a contradiction that we had a lin- 
ear combination of this linear form and all the nontrivial ones that gave (1, 0, . . . , 0) 
as a result. Then we would also have a linear combination of (0, 0, 1, 0, . . . , 0) and 
the rows of the matrix of Figure [9] giving (1, 0, . . . , 0). The coefficients of this linear 
combination would be all zero except from, at most, the ones of the first and the 
third row of the matrix of Figure [9] and the one of (0, 0, 1, 0, . . . , 0). The coefficient 
of the first row would have to be different from zero. But this row has nonzero 
entries in positions were the third row and (0, 0, 1, 0, . . . , 0) have zero entries. Then, 
the linear combination giving (1,0,..., 0) is not possible. 

Also, from the fact that the nontrivial linear forms together with (1,0, . . . ,0) 
and (0, 0, 1, 0, . . . , 0) are linearly independent, we can deduce that there is only 
one possible linear combination of them giving the n-th trivial linear form (the 
coefficients of (1,0,..., 0) and (0, 0, 1, 0, . . . , 0) are and the rest are given by the n- 
th row of the figure of the proof of Lemma[6]) . Then, if we take out of this set of forms 
the first nontrivial one (which had a nonzero coefficient in that linear combination), 
it will not be possible to obtain the n-th trivial form as a linear combination. 

Finally, none of the nontrivial forms can be obtained as a linear combination of 
all but the first and third linear forms, because every nontrivial linear form has a 
nonzero entry either in the first or in the third position. 

These things prove that: 

• Given any trivial linear form different from the n-th one, its complexity is at 
most one because we can divide the rest in the set of trivial and the set of 
nontrivial ones and it will not be in the linear span of any of these sets. 

• Given the n-th trivial form, its complexity is at most one because we can 
construct the next two sets with the rest of the forms: the first one with the 
first and third trivial forms and with all the nontrivial forms except from the 
first one; the second one with the rest of the linear forms. 

• Given any nontrivial linear form, its complexity is at most one because we can 
divide the rest of the forms in the next two sets: the first one with the first 
and third trivial forms and with all the other nontrivial forms; the second one 
with the rest of the trivial forms. □ 
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4 Volumes of polytopes 



We already have systems of linear forms \I/ : Z n ' _2n — > Z™ 2 for every n and 
we know their complexity. Now, we have to choose a suitable convex body K C 
[-N, JV] n2_2n in Conjecture [l] Remember that for our problem we know that Con- 
jecture [l] is a theorem, since we have finite complexity for every n > 3. Recall that 
we are interested in counting the number of re x n Z-magic squares with their entries 
being primes in [0, N]. 

Then, for each n > 3 and N > 0, we will have K = K n {N) = {x e W n2 ' 2n : < 
tpi(x) < N for i = 1, . . . , re 2 }, where the coefficients of ifii are given in figures [2] (for 
n = 3), [5] (for n = 4) and Q (for n > 5), and the term \K n Z" 2 " 2 " n ^(P™ 2 )] in 
Conjecture [T] will be the quantity we are interested in. The first observation is that 
K is convex since it is the intersection of n 2 convex sets. 

We need to calculate the archimedean factor (3^^ = vol n 2„ 2n (i^ D ^ _1 (]R" 2 )). 
Since all the points in our K belong also to ^ _1 (IR™ 2 ), we have /3oo,n = v ol n 2_ 2n (K). 
Since ^ is a system of linear forms, the volume of K n (N) will be equal to jV n2 ~ 2n 
times the volume of the polytope defined by K n (l) = {x G M. d : < ipi(x) < 
1 for i = 1, . . . , n 2 }. So, we only have to compute K n (l) for re > 3. 

Observe that K n (l) is important in its own right. (3oo,n(N) = vol n 2 _ 2n (K n (N)) = 
N n -2n vo i n2 _ 2n (^ n ( 1 )) is t he volume of a convex polytope in M n 2n which points 
with integer coordinates are in bijective correspondence with all the n x re Z- magic 
squares with entries in [0, N]. For large N, the volume of K n (N) will be asymptot- 
ically the same as the number of integer points in it, so the volume of this polytope 
gives an asymptotic for the number ofnxn Z-magic squares with entries in [0, N] 
when iV tends to infinity. We will talk about this in more detail in a moment. 

4.1 Volume for 3x3 Z-magic squares 

We have to calculate the volume of the polytope ^3(1) defined by the next 18 
inequalities (0 < ipi(x) < 1 for i = 1, . . . , 9): < a + b < 1, < a — b — c < 1, 
< o + c < 1, 0<a-b + c<l, < a < 1, 0<a + b-c<l, < a - c < 1, 
0<a + 6 + c<l, < a — b < 1. Each one of these pairs of inequalities defines a 
region bounded by two parallel planes and because we are in a tridimensional space 



we can "see" the polytope, which is the one of Figure 10 
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(1/2,0,1/2) 




(1/2,0,-1/2) 



Figure 10: Polytope K 3 (l). 

Since the polytope is the disjoint union of two pyramids with area of their base 
equal to 1/2 and height 1/2, each pyramid will have volume 1/12 and the polytope 
will have volume 1/6. 

So, /3oo, 3 (A0 = ^ • 

Before jumping to the case n = 4, we could think of a different approach to 
calculate the volume of this polytope, maybe one wich gives a better understanding 
or is more generalizable. An interesting thing to know will be the exact number of 
points with integer coordinates in K = K n {N) for every N . As we were discussing 
before, this will be the exact number of n x n Z-magic squares with entries in [0, N]. 
We have to go to Ehrhart's theory, and look at its fundamental result (see, for 
example, [3] or [5]): 

Theorem 8. (Ehrhart's Theorem) If P is a convex polytope whose vertices have 
rational coordinates, then the number of points with integer coordinates in the dila- 
tions NP for N = 0, 1, 2 . . . is a quasi-polynomial in N whose degree is the dimension 
of P and whose period divides the least common multiple of the denominators of the 
vertices of P. 

In this case, since the least common multiple of the denominators of the vertices 
of ^3(1) is 2, we will have that the Ehrhart quasipolymial has degree 3 and period at 
most 2. A direct way to obtain the quasipolynomial, Es(N), would be to interpolate 
it. For that we would need (3 + 1)2 = 8 values. E^(N) is the number of points 
with integer coordinates, x, satisfying < tpi(x) < N for i = 1, . . . , 9 and, with a 
computer (or even by hand!), it is easy to compute: 

£ 3 (0) = 1 £ 3 (1) = 2 £ 3 (2) = 7 £ 3 (3) = 12 
£ 3 (4) = 25 £3(5) = 38 £3(6) = 63 £ 3 (7) = 88 
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Now, we can interpolate our quasipolynomial 

f iV 3 /6 + N 2 /2 + AN/3 + 1 if N = (mod 2) 

E-i(N) = \ 

{ N 3 /6 + N 2 /2 + 5N/6 + 1/2 if N =1 (mod 2) 

which gives the exact number of points with integer coordinates in K^(N), i. e. the 
exact number of 3 x 3 Z- magic squares with entries in [0, N]. The coefficient of the 
leading term is 1/6, exactly the volume of K^(l). This is not a coincidence, and 
it is a general fact that the leading coefficient of the Ehrhart quasipolynomial of a 
polytope is its volume. It is an inmediate consequence of the observation that the 
volume of a d-dimensional polytope can be thought as the limit when N — > 00 of 
the number of integer points in its iV-dilation divided by N d . With all this in mind, 
we can go to the next case. 

4.2 Volume for 4x4 Z-magic squares 

Our approach is now clear. We want to calculate the volume of the polytope 
K,4(l) given by the 16 pairs of inequalities given by < tj^i(x) < 1 for i = 1, . . . , 16, 
where the coefficients of ipi(x) are given by the z-th row of the matrix of Figure |8j 
We know that the volume is given by the leading coefficient of the quasipolynomial, 
E±(N), that gives the exact number of integer points satisfying < ipi(x) < ./V for 
i = 1,...,16. 

We will compute E^(N) by interpolation. The period of E^N) divides the least 
common multiple of the denominators of the vertices of K^(N). So, in order to 
know how many values of E^N) we need to interpolate it, we have to calculate 
the vertices of the polytope. There are different methods to compute the vertices 
of a polytope (see, for example, |13| or [1] for surveys of vertex finding algorithms). 
The most direct of them would be to look at all the intersections of 8 of the 32 
equalities ipi(x) = or ipi(x) = 1 for % = 1, . . . , 16 consisting of a single point and 
then checking if the point satisfies all the inequalities < tpi(x) < 1 for i = 1, . . . , 16. 
You do not want to do this by hand, but your computer will tell you that the 178 
vertices of -^4(1) are the ones listed in Appendix [Aj The only important thing for 
us is that the denominators of the coordinates of the vertices are one, two or three, 
so their least common multiple is 6. Then we are looking for a quasipolynomial, 
E4(N), of degree 8 and period dividing 6. In order to interpolate E^N) we would 
then need (8 + 1)6 = 54 values. 

In Appendix [B] we calculate the 54 needed values. With them we interpolate our 
quasipolynomial and we have 
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' 83897V 8 , 83897V 7 , 55317V" , 108777V 5 , 86637V 4 , 14371/V 3 , 1437V 2 , 10677V , i 
120960 i " 15120 i " 2592 2160 + 1080 + 1620 + 21 + 315 +1 

if N = (mod 6) 



83897V 8 , 83897V 7 
120960 i " 15120 

if JVsl (mod 



55317V 6 , 108 777V r ' , 691697V 4 
2160 



570797V 3 , 43037V 2 



2592 



8610 



6 180 



672 



+ - 



5010 



128 



83897V 8 
120960 



83897V 7 
15120 



55317V" i 108 777V 5 
2592 "r 2160 



86637V 4 , 143717V 3 , 1437V 2 
1080 + 1620 + 21 



10677V I 97 
315 81 



if N = 2 (mod 



S 4 (iV) = < 



83897V 8 
120960 



83897V' 
15120 



55317V 1 ' 
2592 



108777V ■' 
2160 



691697V 4 
8640 



570797V 3 
6480 



43037V 2 
672 



136077V 
5040 



37 

128 



if N = 3 (mod 6) 



83897V 8 i 83897V 7 , 
120960 " 1 " 15120 ~* 

if N = 4 (mod f 



55317V" i 108777V 5 , 86637V 4 , 143717V 3 , 1437V 2 , 10677V , -, 
2592 + 2160 + 1080 + 1620 + 21 + 315 + 



83897V 8 
120960 



if 



N = 5 



83897V 7 , 55317V" , 108 777V 5 , 691697V 4 , 570797V 3 , 43037V 2 , 136 077V , 5045 
15120 "r" 2592 "r 2160 "f 8640 " + " 6480 " + " 672 "r 5040 " + " 10368 

(mod 6) 



which gives the exact number of points with integer coordinates in K^N), i. e. 
the exact number of 4 x 4 Z-magic squares with entries in [0,7V]. The coefficient 
of the leading term, 8389/120960, is the exact volume of K±(l) (one can confirm 
this value, for example, with polymake; see |http : //www ."p olymake . org/[). You may 
wish to compare the quasipolynomials E^(N) and E^{N) with the quasipolynomials 
of period 2 computed in [lj, depending on the magic sum and not on a bound for 
the entries. Then, ^(N) = Sj. 



4.3 Volume for higher values of n 

It is reasonable to think that if for the case 3 x 3 we were able to interpolate our 
quasipolynomial by hand and for the case 4 x 4 we could do it with some computer 
time, then for the case 5 x 5 we would be able to interpolate the corresponding 
quasipolynomial with just a little more computer time. 

But this is not exactly the case. For the case of 5 x 5 Z-magic squares, we would 
have a polytope, K^{1), in a 15-dimensional space defined by 25 pairs of inequalities 
given by < ipi(x) < 1 for i = 1, . . . , 25. It is easy and fast to find some vertices of 
this polytope with denominators 3, 5, 7 and 8, for example 

(0,0,0,1,1,0,1,1,0,1,0,1,0,1,1) (0,0,0,l,l,l,0,l,0,0,|,f,0,l,|) 
(0,0,0,l,l,0,f,l,f,f,l,f,0,l,f) (0,0,1,0,1,1,|,0,|,1,0,|,|,0,|). 

This means that the least common multiple of the denominators of the vertices of 
K§(1) is at least 3 • 5 • 7 • 8 = 840. To interpolate a quasipolynomial of degree 15 and 
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period dividing 840 we would need (15 + 1)840 = 13440 values. In order to calculate, 
for example, the largest one of these values, in principle, one would have to consider 
the points with integer coordinates in [0,13339] 15 and check if the 20 inequalities 
given by the nontrivial linear forms are satisfied. Even considering that we could 
reduce the number of values to roughly half (6719 values) with Ehrhart-Macdonald 
Reciprocity Law (see Appendix JbJ) , and that with efficient programing the number 
of calculations can be further reduced, one would have to make a very very high 
number of checks. And of course this grows extremely quickly with n. 

It is true that we could try finding the quasipolynomial with other ideas and 
also that for our problem we do not need the quasipolynomial E n (N) but only its 
leading coefficient, which gives the volume of K n (l). And it is also true that there 
are many other ways to calculate the volume of a polytope. But, to the best of 
our knowledge, no known method would give a formula for the exact volume of 
K n (l) for general n, and giving the exact value of the volume would be hard even 
for relatively small values of n. As an example illustrating this, for the Birkhoff 
polytope, one of the most important and studied polytopes, which is a similar but 
definitely simple^ polytope than the one we are considering, the volume is only 
known for n < 10 and the case n = 10 took almost 17 years of computer time (see 
[2] and http : / / www . math . binghamt on. edu/dennis/Birkhof f /[ ) . 

So we cannot give a formula for the exact volume of the polytope K n (l), but 
we can prove that the volume will be a nonzero rational number. The fact that it 
will be a rational number is easy to see with our approach because the coefficients 
of the polynomials we are interpolating are given by the solutions of linear systems 
with integer coefficients. 

To show that the volume is never zero, we first make the next observation: 

Lemma 9. Consider the linear forms given by the elephant basis for some n > 5. 
Then, for every i = 1, . . . , n 2 the sum of the coefficients of the linear form ipi is 1. 

Proof. Since the sum of magic squares gives a magic square, when we add up the 
n 2 — 2n vectors of our elephant basis we obtain a magic square. The obtained magic 
square has 1 as an entry in all the positions of the skeleton of the elephant basis. 

The Birkhoff polytope B n is the convex polytope in M whose points are the doubly stochastic 
matrices, i.e. the n x n matrices whose entries are non-negative real numbers and whose rows and 
columns each add up to 1. B n is a convex polytope with integer vertices (which are well known). 
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Then, this magic square is the only magic square with ones in all that positions, so 
it has to be the magic square with n 2 ones as entries. Since the coefficients of ijji are 
the entries on the i-th position of the vectors of the elephant basis, we are done. □ 

Then, the point (A, . . . , |) G M n2 ~ 2n belongs to K n (N) because ^(|, . . . , \) = 
1/2 for every i = 1, . . . ,n 2 . Since all the coefficients of all our linear forms have 
absolute value less than or equal to n (see, for example, Figure [7]) then all the points 
(I, ■ • • > I) + w e j for j = 1, ... n 2 — 2n, where Cj is the j-th vector of the canonical 
basis of the (n 2 — 2n)-dimensional Euclidean space, belong to K n (l). Because the 
vectors joining (A, . . . , A) to these n 2 — 2n points are linearly independent, we have 
that the convex hull of the n 2 — 2n + 1 points has nonzero volume, and because 
K n {N) is convex we have that it contains that convex hull. Then, K n {l) has nonzero 
volume. 

Since the volume of K n {\) is less than or equal to one (in fact it is much smaller 
than one) because the inequalities given by the trivial linear forms define the unit 
cube, we have proved: 

Lemma 10. For every n > 5, /3 00jn (A r ) = c{n)N n2 ~ 2n where c(n) G (0, 1] is the 
rational number that gives the volume of K n {l). 

5 The local factors 

The only ingredients we are missing in order to give our asymptotics are the 
local factors f3 P)n . Recall that for each prime p, 

9 

rr 
i=l 

where is the local von Mangoldt function, A% p : Z p — > M + , defined by A^ p (0) = 
and A Zp (6) = ^ y = ^ T if6/0. 

Then, the important point is determining for which elements m £ Z™ 2 ~ 2n , we 
have ipi(m) ^ (mod p) for every i = 1, 2, . . . , n 2 . In order to do that we will use 
the inclusion-exclusion principle. As we will observe, for small primes we can have 
a different behaviour but one can always determine a nice formula giving the value 
of /3 p for every p from some point on. 
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5.1 Local factors for n = 3 

We can think we have the system of 9 linear forms \I/ : H? p — > defined by 
^(a, b, c) = (a + 6, a — b — c, a + c, a — 6 + c, a, a + 6 — c, a — c, a + 6 + c, a — 6), and we 
want to determine how many of the p 3 elements of Zp 1 give nonzero values modulo 
p for all nine linear forms. 

We use the inclusion-exclusion principle. There are p 3 elements in Zp 1 . If none of 
the linear forms is trivial in Z^, any one of them is cancelled for p 2 values of Zp\ We 
have to be careful because even if none of our linear forms is trivial in the integers, 
it could be the case that some of them were trivial when reducing modulo p for some 
small values of p. Since all the coefficients of our linear forms are zeros, and plus or 
minus ones, this is not the case. 

Now, if none of the linear forms is a multiple of another, any of the (2) systems 
formed by two of the linear forms is cancelled for p values of Zp\ This is satisfied 
for every p > 3 but for p = 2 some of the linear forms are identical to others (for 
example, the second one is the same as the fourth, the sixth and the eighth). We 
will later study the case p = 2 separately. 

The next step is computing, for every system of three of the linear forms, the 
rank of the matrix formed with their coefficients as rows. Again, small primes can 
be tricky. For example, if we select the last three linear forms, the corresponding 
matrix with integer coefficients and its Hermite normal fornF] are respectively: 



/ 1 -1 ^ 


\ 


f 1 





2 \ 


1 1 1 


and 





1 


2 


\ 1 -1 ) 




\o 





3 / 



Then, the matrix in our example has rank 2 in Z3 and has rank 3 in Z p for every 
p > 3. We will also later study the case p = 3 separately. It makes sense to compute 
the Hermite normal form of all the matrices with the coefficients of 3, 4, 5, 6, 7, 8 or 

3 There are different conventions to define the Hermite normal form of a matrix with coefficients in 
Z. For example, we can allow the next three unimodular elementary row operations: interchanging 
two rows, multiplying one row by —1 and adding an integer multiple of a row to another one. Then, 
the Hermite normal form of a matrix with coefficients in Z can be defined as the unique matrix 
obtained from it by unimodular elementary row operations such that the first nonzero entry of each 
row is positive and strictly to the right of the first nonzero entry of the row on top of it and such 
that all the entries in the same column than the first nonzero entry of each row are less than it and 
greater than or equal to 0. 
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the 9 forms as rows. If we then look at the largest of all the first nonzero elements 
of a row we will see that it is equal to 4. Then, we know that for every prime p > 5 
the Hermite normal form of these matrices in Z has exactly the same rank when we 
reduce its coefficients modulo p. This means that we can study at the same time 
/9 P) 3 for every prime p > 5. 

We have that 8 of the matrices with three of the linear forms as rows have rank 
2 and the other 76 have rank 3 in Z. Also, any matrix with 4, 5, 6, 7, 8 or the 9 
forms as rows has rank 3 in Z. This means that we can use the inclusion-exclusion 
principle to deduce that for every prime p > 5 the number of elements m S Z 3 , with 
tpi(m) (mod p) for every i = 1, 2, . . . , 9 is: 

p 3 - 9p 2 + 36p - (8p + 76) + 126 - 126 + 84 - 36 + 9 - 1 = p 3 - 9p 2 + 28p - 20. 

Then, for every prime p > 5 

_ p 3 - 9p 2 + 28p - 20 / p \ 9 
p ' p 3 \p — 1 J 

For p = 2, the four elements of Z 3 starting with cancel the fifth linear form. 
The element (1,0,0) does not cancel any of the nine linear forms. The other three 
elements starting with one either cancel the first or the third linear form. Then, 
02,3 = gs • 2 9 = 2 6 . Similarly, for p = 3 it is easy to see that the only two elements of 

n 

P3,3 - p • [2) -38- 



3 which do not cancel any of the nine linear forms are (1, 0, 0) and (2, 0, 0). Then 

06 



5.2 Local factors for n = 4 

In this case we have the 16 linear forms which coefficients are given by the rows 
of the matrix of Figure [8j We will determine how many of the p s elements of Z^ 
give nonzero values modulo p for all the sixteen linear forms. 

First of all we compute the Hermite normal form of all the matrices formed with 
some of the rows of the matrix of Figure |8j The largest of all the first nonzero 
elements of a row appearing in all these matrices is 3, so, as before, for every prime 
p > 5 the Hermite normal form of these matrices in Z has exactly the same rank 
when we reduce its coefficients modulo p. 

Now we compute the ranks in Z of all these matrices and we have: all 16 one 
row matrices have rank 1, all 120 two row matrices have rank 2 and all 560 three 
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row matrices have rank 3; 10 four row matrices have rank 3 and 1810 have rank 4; 
120 five row matrices have rank 4 and 4248 have rank 5; 708 six row matrices have 
rank 5 and 7300 have rank 6; 32 seven row matrices have rank 5, 2656 have rank 

6 and 8752 have rank 7; 433 eight row matrices have rank 6, 6553 have rank 7 and 
5884 have rank 8; 32 nine row matrices have rank 6, 2656 have rank 7 and 8752 
have rank 8; 708 ten row matrices have rank 7 and 7300 have rank 8; 120 eleven 
row matrices have rank 7 and 4248 have rank 8; 10 twelve row matrices have rank 

7 and 1810 have rank 8; all 560 thirteen row matrices have rank 8 and the same 
thing happens with the 120 fourteen row matrices, the 16 fifteen row matrices and 
the only sixteen row matrix. By the inclusion-exclusion principle, for every p > 5 
the number of elements m S Z 8 , with tpi(m) ^ (mod p) for every i = 1,2, ... ,16 
is: 

p 8 - 16p 7 + 120p 6 - 550p 5 + 1690p 4 - 3572p 3 + 5045p 2 - 4257p + 1539. 
Then, for every prime p > 5 

p 8 - 16p 7 + 120p 6 - 550p 5 + 1690p 4 - 3572p 3 + 5045p 2 - 4257p + 1539 / - x 16 



A 



p.i 



p-1 



There is 1 element of Z| and there are 34 elements of Z| not cancelling any of 
the sixteen linear forms, so 

1 16 s 34 /3\ 16 17-3 8 

02,4 = ^ • 2 16 = 2 8 and fo, 



2 8 " — ^ 38 ^2/ 2 15 

5.3 Local factors for higher values of n 

In the general case, it is clear that we can follow the same steps. Since every 
two of our linear forms are linearly independent, for primes p > po(n) we will have: 

p" 2 -^ - n^ 2 - 2 "- 1 + P< n2 _ 2n _ 2 (p) { p 

Pi 



P,n pn 2 -2n \p _ I 

2n 2 -2n 2 2n 2 -2n-l , p / \ 

P — U p + f<2n 2 -2n-2\P) 

r .2n 2 —2n ri 2 ri 2n 2 — 2n—l i p „ f n \ 
P —np + ^=2n 2 -2n-2\P) 

1 P <2n 2 -2n-2(p) 
P=2n 2 -2n{P) 

where P<d{p) and P=d(p) are (possibly different in each occasion) polynomials in 
p of degree less than or equal to D and D respectively. 

The existence of the logarithm function makes the theory of infinite products 
essentially equivalent to the theory of infinite series. In particular, it is well known 
that (see, for example, [TT]): 
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Lemma 11. If {a,}^ 1 is a sequence of real numbers such that < |aj| < 1 and 
YliLi \ a i\ < oo then n£i(l + a converges to a nonzero real number. 

We know that for every prime, p, we have /3 p>n > because m = (1, . . . , 1) G 
Z n2_2n makes ipi(m) = 1 for every i = 1,2,... ,ra 2 . Also, since there is a constant, 
C(n) > 0, such that for p > po(n) we have 

P<2n 2 -2n-2{P) 
P=2n 2 -2niP) 

then, by the previous lemma, we deduce: 

Lemma 12. For every n > 3, Y\ p prime (3 p , n converges to a strictly positive real 
number. 



< 



C(n) 



6 Asymptotics for Magic Squares of Primes 

The work of the previous sections and the results mentioned in the Introduction 
prove the next three theorems: 

Theorem 13. The number of 3 x 3 Z-magic squares with their entries being prime 
numbers in [0, N] is 

N 3 

(l + 1 )6 3 -^-, 
log N 

where 

8 A A p' 3 \P — 1/ 

p prime x 7 

p>5 

Observe that, apart from the asymptotics for 5-term arithmetic progressions of 
primes, this is one of the first "natural" applications of the work of Green, Tao and 
Ziegler to a system of linear equations of complexity 3 (see Section [3]). 

For n = 4 we have complexity 1 (see Section [3]) and the asymptotic: 

Theorem 14. The number of 4 x 4 7L-magic squares with their entries being prime 
numbers in [0, N] is 

N 8 
log 16 iV' 
where 

16 



(1 • "(I))©-,, 



p~ _ 34654959 , , 
° 4 ~ 573440 IIP p™ne 



p>5 

76.758. 



n p 8 -16p 7 + 120p 6 -550p 5 + 1690p 4 -3572p 3 +5045p 2 -4257p+1539 / p \ 
P prime p s \p—l J 
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^yn 2 — 2n 



Finally, for n > 5 (see Section [3]) we have complexity 1 and the "order of mag- 
nitude of the asymptotics" : 

Theorem 15. The number ofnxn Z-magic squares with their entries being prime 
numbers in [0, N] is 

(l + o(l))6 r 

log N 

where & n = vol n 2_ 2n {K n (l)) ]J p prime l3p,n is a nonzero real number. 

7 Magic squares with different entries 

Up to this point we have not been careful with the fact that there could be 
repetitions among the entries of our Z-magic squares. In this section we prove that 
the asymptotics for Z-magic squares of primes with different entries are exactly the 
same as the asymptotics for Z-magic squares of primes. We first observe: 

Lemma 16. For every n > 3, there are well known constructions ofnxn Z-magic 
squares with their n 2 entries being the first n 2 natural numbers (these are usually 
called normal magic squares). 

Maybe the best known constructions for normal magic squares of odd side are 
the ones popularly known today by the names of Simon de la Loubere and Bachet 
de Meziriac. For normal magic squares of even side there are also different construc- 
tions, for example those of Devedec to name some. All these constructions can be 
found in [12] and many others very easily searching on the Internet. 

The only important thing for us is that there exist magic squares with all their 
entries being different. With this in mind we can prove the next lemma. 

Lemma 17. For every n > 3 the number ofnxn Z-magic squares with their entries 
being primes in [0, N] is asymptotically equal to the number ofnxn Z-magic squares 
with their entries being different primes in [0, N] when N — > oo. 

Proof. In the same way that we proved that the 2n equations defining nx n magic 
squares were linearly independent exhibiting magic squares that satisfied some of 
the equations but not the others (see pa ge [5] ), the fact that there are Z-magic squares 



with all their entries different (Lemma 16 ) proves that any of the ) equations of 



the type xi — Xj = for i £ [1, n 2 — 1] and j £ [i + 1, n 2 ] is linearly independent with 
the In equations defining nx n magic squares (see page [4]) . 
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Then, the number of Z-magic squares with their entries being different elements 
of [0,N] is O^iV" 2 - 2 "- 1 ). This is certainly o n (N n2 ~ 2n /log 11 ' 2 N) and then we are 
done. □ 
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Vertices of K^l) 

The 178 vertices of ^4(1) in lexicographic order are: 



0,0,0,0,0,0,0,0) 

i 1 i i 1 11 

u , u 5 2 ' ' 2 ' 2 ' 2 ' / 

1 J - 1 1) 

0,0,1,1,0,1,0,1) 
0,0,1,1,1,1,0,0) 
0,1,1,0,0,|,0,|) 
i 1 1 1 i I 1) 

0,1,0, |,1,0, |,0) 
0,1,0,1,0,1,0,1) 

0,1,0,1,1,1,1,1) 

Oillnllll 

u , x > 2 ' 2 ' u > 2 ' 2 > / 

0,1,1,0,0,1,1,1) 
1 1 1 1 1 1 I) 

|, 0,1, o,o,|,o, |) 
1 
3 
1 
2 
1 
2 
1 
2 
1 
■> 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
1 
2 
2 

:<> 



1 2 2 x : 2 2n 

J-j 3; 3, J-j 3; 3) 

0,0,1,1,1,0,0) 

o,|,|4,i,o,o) 

0,1,0,0,|,0,0) 

1 I 1 - -) 
0,1,|,1,0,0,0) 
0,1,1,1,1,0,1) 

io,|,o,o,i,i) 

±,0,1,1, ±,1,0) 
ill 10ii) 

2' ' 2''' 2' 2 ^ 

1,0,0,0,1,1,0) 

1 1 1 1 M 

1,0,1,1,0,0,0) 
1,0,1,1,|,1,|) 
1 I I 1 1 1) 

' 2 ' 2 ' ' 2 ' ' / 

1,1,1,1,1,1,1) 
1 1 I 1 1 I I) 

x , ^, 2 ' ' ' 2 ' 2 / 

? 1 I ? 1 1) 

3 ' 3 ' ' 3 ' > / 

1,0,0,0,0,0,1,0) 
1,0,0,1,0,1,1,1) 

1 1 i 1 1) 



1 ' ^ ' 2 2 



1,0,1,0,0,0,1,1) 
1,0,1,0,±,0,±,0) 
1,0,1,1,0,1,1,1) 
1,0,1,1,1,1,1,1) 



0,0,0,1,0,1,0,1) 
0,0,1,0,0,±,0,±) 

0,0,14,0,1,0,1) 

0,0,1,1,|,1,|,1) 
0,0,1,1,1,1,0,1) 
0,|,1,0,0,|, |,1) 
0,1,0,0,0,0,0,0) 

0,1,04,1,0,0,0) 
0,1,0,14,14,1) 

0,1,0,1,1,1,0,0) 

1 i I 1 I I 0] 

w , 1 i 2 ' 2 ' ' 2 ' 2 ' u / 

0,1,1,0,1,0,0,0) 
0,1,1,1,1,1,0,1) 

14,0,1,14,0,0) 

1,0,0,1,0,044) 

0, §,0,0,0,0,0) 

044,14,0,0) 
0,1,0,04,04) 
0,14,0,1,1,1) 

1 I 1 - -) 

, , 2' ' ' 2 ' 2/ 

0,1,1,144,1) 

± 1 1 1 ±) 

2' 2 ' ' ' 2 ' 2 / 

1 1 - - 1) 

1,14,1,1,0,0) 
1,0,0,04,14) 
1,04,0,1,1,1) 

1 1 1 1 -) 

J -i vy 7 2' ' ' 2' 2 

1,0,1,14,1,1) 

1 1 1 I 1 1) 

1,1,0,04,1,1) 
± - ± ±) 
i,o,i,i,|,i, |) 
1,0,04,0,044) 

1,0,0,1,1,0,0,0) 

1 n I 1 1 I i n] 

x , u , 2 , 2 ' ' 2 ' 2 ' / 

1,0,1,0,0,14,0) 
1,0,1,0,1,0,1,0) 

1,0,1, |,o,i, 
14,0,0,044,0) 



0,0,0,1,1,0,0,0) 
0,0,1,0,0,1,0,0) 

0,0,14,0,14,1) 
0,0,1,1,14,04) 
04,04,1,0,0,0) 
04,14,0,1,0,1) 

0,1,0,0,0,0,1,1) 

0,1,04,1,04,0) 

0,1,0,1,1,0,1,1) 

1 i 1 1 1 0) 

u , 2 ' ' 2 ' 2 ' 2 ' / 

0,14,1,1,1,1,1) 
0,1,1,0,144,0) 

0,1,1,1,1,1,1,1) 

12i 2 2 1 2 1 \ 

3 ' 3 ' > 3 ' 3 ' 3 ' 3 ' / 

1 n ± 1 i 01 

2 , u, u, 2 , 2 , 2 , u , u / 

1 n I n) 

2 ' ' 2 ' 2 ' ' 2 ' ' / 

1,0,1,1,0,1,0,1) 
|,0,1,0, 1,0,0,0) 

1,0,1444,0,0) 

1,0,14,1,1,0,0) 

- - - i 01 

2 , 2 , u , u , w , 2 ' 2 ' / 

14,04,1,0,0,0) 

n 1 1 1 1 0) 

2 , 2 ' ' ' ' 2 ' 2 ' ' 

1 I 1 I 1 1 I i) 

2 , 2 ' > 2 ' ' > 2 ' 2 / 

|,1,0,0, 1,0,1,0) 

1 1 i i i 01 

2 ' ' u ' 2 ' 2 ' 2 ' u ' u > 

1,1,04,1,0,1,0) 
1,14,0,0,0,1,1) 

1 1 1 a 1 1 11 

2 ' ' 2' 2' 2 ' ' ' / 

1,1,1,0,1,14,0) 



-011-- - 1) 
3 ' ' ' ' 3 ' 3 ' 3 ' ' 



2 1 21 1 21qN 

3 ' > 3 ' 3 ' ' 3 ' 3 ' / 

1 1 1 1 1 01 

, , , 2 ' 2 ' 2 ' 2 ' / 
1,0,0,1,144,0) 

1,04,1,0,1,1,1) 
1,0,1,0,04,14) 
1,0,14,0,144) 
1,0,1,14,14,1) 
14,04,0,0,1,1) 



a n I I 1) 

u , u , 2 ' 2 ' ' 2 ' 2 ' / 

0,0,1,0,1,0,0,0) 

0,0,14444,1) 

0,0,1,1,144,1) 

04,0,1,14,04) 

04,14,1,1,0,0) 

0,1,0,04,04,0) 

0,1,04,1,044) 

0,1,0,1,14,04) 

0,14,0,1,04,0) 

0,1,1,0,044,1) 

l l U I 1 1) 

, , , 2' 2' 2' 2' ' 

1 n 1 2 q 1 2 



0, 



3 , 3 , "j 3 , 3 , 



1) 



1,0,0,14,1,0) 

0,0,1,0,14,1) 

0,144,0,0,0) 

0,1,1,1,1,0,0) 

0,1,1,0,1,0,1) 

1 i - i J J) 

u , 1 i 2 ' 2 ' 2 ' ' / 
0,1,14,1,0,1) 

- 1 - -) 

2 , v, 2 , u , u , 2 ' 2 / 
1,0,1,0,14,1) 

1,14,0,1,1,1) 

1 1 1 1 1 1 1) 

2 , a j x , 2 ' 2 ' / 

1,0, |,0,0,1,1) 
1 1 1 - 1 1) 

, , 2 ' 2 ' 2 ' ' / 
1,0,14,1,1,1) 
1,1,0,1,0,1,0) 
I44,l,i0.0) 



1,1,1, 

,0, 



,0) 



3,3,", 3,3,3,3 
1,0,0,0,0,0,0,0) 

1,0,0,1,044,1) 



1,04,0,0,04,0) 

1 1 1 i 1 1 1) 

2 ' ' 2 ' 2 ' 2 ' / 

1,0,1,0,0,1,0,0) 

1,0,14,0,14,1) 

1,0,1,1,1,1,0,0) 

14,04,1,0,1,0) 
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(1, |,0,1 5 1, |, |,0) (1,1,0,1,1,1,1,1) (1,1,1,0,0,1,1,1) (1,1,1,1,0,1,1,1) 

(1,1,0,0,0,0,1,0) (1,1,0,0,0,0,1,1) (1,1,0,0,0,1,1,0) (1,1,0,0,0,1,1,1) 

(1,1,0,0,1,0,1,0) (1,1,0,0,1,0,1,0) (1,1,0,1,1,1,1,0) (1,1,0,1,1,0,1,0) 

(1,1,0,1,1,0,1,0) (l,l,0,|,l,i,|,0) (1,1,0,1,0,1,1,1) (1,1,0,1,1,0,1,1) 

(1,1,0,1,1,1,1,1) (1,1,1,0,1,1,1,0) (1,1,11,1,1,1,0) (1,1,1,0,0,1,1,1) 

(1,1,1,0,1,0,1,0) (1,1,1,1,1,1,1,1) 



These vertices were obtained with the algorithm described in Section 4.2 Ex- 
actly the same vertices are obtained using cddlib, the implementation of the double 
description method of Motzkin, Raiffa, Thompson and Thrall (Copyright by Komei 



Fukuda, http : //www, if or .math, ethz. ch/~fukuda/cdd_home/ cdd.htmlj ). 



B 54 values of E 4 (N) 

Recall that E^(N) is a quasipolynomial giving the exact number of integer points 
that satisfy < ipi{x) < N for i = 1, . . . , 16, where the coefficients of tpi(x) are given 
by the i-th row of the matrix of Figure [8j Since the degree of E^(N) is 8 and its 
period divides 6, we need 54 values of E^{N) in order to interpolate it. 

Computing each one of these values essentially consists of checking for each one 
of the integer points in [0, N] 8 if it satisfies the 16 inequalities given by the eight 
nontrivial linear forms being between and N, including both of them. It is true 
that we can be clever and reduce considerably the number of checkings, but this 
gives us an idea of the order of magnitude of the number of checkings needed. With 
this in mind, it would be nice if we could reduce the number of values of E^N) 
needed to interpolate it. The next fundamental result (see, for example, [3] or |14j ) 
allow us to reduce significantly that number of values, from 54 to roughly the first 
half of them. 

Theorem 18. (Ehrhart-Macdonald Reciprocity Law) Suppose that P is a 
convex polytope whose vertices have rational coordinates. Then the evaluation of its 
Ehrhart quasipolynomial, E{N), at negative integers yields 

E(-N) = {-l) dimP E°{N), 

where E°(N) is the number of integer points in the interior of the dilation NP. 

How to obtain the values of E° from the values of El E°(N) is the number 
of points with integer coordinates in the interior of the polytope K n (N) = {x £ 
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R n ~ 2n : < tpi(x) < N for i = 1, . . . ,n 2 }. A point is in the interior of K n {N) if 
and only if it satisfies those 2n 2 inequalities but none of the corresponding equalities, 
so K°(N) = {x £ R n2 - 2n : < ipi(x) < N for i = 1, . . . , n 2 } = {x G IR™ 2-2 ™ : 1 < 
M x ) < N- 1 for i = 1,... ,n 2 }. 

Because of the property that for every i = 1, . . . , n 2 the coefficients of ^ add up 
to 1 (see Lemma [9]), substracting 1 to all of the coordinates of the elements of K°(N) 



gives a biyection between this set and {x G M™ 2 2n : < ViO^) < -/V — 2 for i 



1, . . . , n 2 } = K n (N - 2) for iV > 2. Then, Theorem 18 implies that m our case 



E 4 (-N) = E 4 (N - 2) for N > 2. 



This, and a little program checking the number of points x G [0, N] 8 that satisfy 
< ipi(x) < iV for the 2n nontrivial linear forms and iV G [0, 26], together with the 
fact that E°(l) = 0, gives: 



E 4 (-l) 


= 


E 4 (0) = 


= 1 = E 4 {-2) 


E 4 (l) = 


34 = E 4 (-3) 


E 4 {2) = 


-- 621 = E 4 (-A) 


E 4 (3) = 


5400 = E 4 (-5) 


£4(4) = 


= 30277 = E 4 (-6) 


E 4 (5) = 


125794 = E 4 {-7) 


£4(6) = 


= 423097 = E 4 (-8) 


£4(7) = 


1214992 = E 4 {-9) 


£4(8) = 


= 3089369 = E 4 (-10) 


£4(9) = 


7130034 = £4 (-11) 


£4(10) 


= 15210869 = E 4 (-12) 


£4(11) = 


= 30399592 = E 4 (-13) 


E 4 {12) 


= 57508653 = E 4 (-U) 


E 4 (13) = 


= 103807042 = #4 (-15) 


E 4 {U) 


= 179946753 = E 4 (-16) 


E A (1S) ~- 


= 301109616 = E 4 (-17) 


£ 4 (16) 


= 488451089 = £ 4 (-18) 


E 4 (17) -- 


= 770830866 = E A (-19) 


£ 4 (18) 


= 1186938765 = E 4 (-20) 


E 4 (19) -- 


= 1787779544 = E 4 (-21) 


£4(20) 


= 2639668773 = E 4 (-22) 


£4(21) = 


= 3827663858 = £4 (-23) 


E 4 (22) 


= 5459641001 = E 4 (-24) 


£4(23) = 


= 7670885920 = E 4 (-25) 


£ 4 (24) 


= 10629486297 = E 4 (-26) 


^4(25) = 


= 14542317074 = E 4 (-27) 


£4 (26) 


= 19662006197 



Although in this case we can directly calculate with a computer the 54 values 
E 4 (0), ... , £"53(0) -and we have done so!-, this takes around 217 times the time it 
takes to calculate the above values. This explicitly shows the "power" of Ehrhart- 
Macdonald Reciprocity Law. 
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